Multidimensional solitons in periodic potentials 
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Abstract 

The existence of stable solitons in two- and three-dimensional (2D and 3D) 
media governed by the self-focusing cubic nonlinear Schrodinger equation with 
a periodic potential is demonstrated by means of the variational approxima- 
tion (VA) and in direct simulations. The potential stabilizes the solitons 
against collapse. Direct physical realizations are a Bose-Einstein condensate 
(BEC) trapped in an optical lattice, and a light beam in a bulk Kerr medium 
of a photonic-crystal type. In the 2D case, the creation of the soliton in a 
weak lattice potential is possible if the norm of the field (number of atoms in 
BEC, or optical power in the Kerr medium) exceeds a threshold value (which 
is smaller than the critical norm leading to collapse). Both "single-cell" and 
"multi-cell" solitons are found, which occupy, respectively, one or several cells 
of the periodic potential, depending on the soliton's norm. Solitons of the 
former type and their stability are well predicted by VA. Stable 2D vortex 
solitons are found too. 
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Typeset using 



Recent successful observation [1,2] and theoretical consideration [3,4] of one-dimensional 
(ID) solitons in self- attractive Bose-Einstein condensates (BECs) with a negative scattering 
length raises a question whether higher- dimensional solitons can be created in BECs [5]. 
One possibility, recently proposed in Refs. [6,7], is to induce self-trapping of a 2D conden- 
sate by means of ac external magnetic field via the Feshbach resonance, which makes the 
scattering length a periodically sign-changing one (a similar mechanism was earlier shown 
to stabilize light beams in a bulk optical medium consisting of periodically alternating self- 
focusing and defocusing layers [8]). On the other hand, an efficient technique to control 
BEC dynamics is based on optical lattices (OLs), i.e., interference patterns generated by 
laser beams illuminating the condensate (see earlier works [9] and recent ones [10]). A possi- 
bility to create higher-dimensional gap solitons in self- repulsive BECs loaded into 2D or 3D 
OL was proposed in Ref. [11], where it was shown that soliton-like structures emerge due to 
modulational instability of BEC in the periodic potential (a similar mechanism for the ID 
case was investigated in Refs. [12,13]). Nonlinear localization of atomic Bloch waves in the 
form of 2D matter-wave gap solitons was shown by numerical solution of the corresponding 
eigenvalue problem [14]. 

The objective of the present work is to show that stable 2D and 3D solitons can be 
created in a more straightforward (and more relevant to the experiment) way, namely, in 
self-attractive BECs trapped in OLs. In the absence of the external potential, stationary 
soliton solutions to the corresponding self-focusing nonlinear Schrodinger (NLS) equation 
can be easily found, but they are unstable due to the, respectively, weak and strong collapse 
in the 2D and 3D NLS equation [15]. We demonstrate that periodic potentials of the OL 
type can stabilize 2D and 3D solitons of various types (in particular, both fundamental and 
vortex solitons are found to be stable in the 2D case). 

The same model finds another physical realization in the 2D case, in terms of a light beam 
propagating in a bulk self-focusing Kerr medium with transverse modulation of the refractive 
index, i.e., a medium of a photonic-crystal (PhC) type (see, e.g., Ref. [16]). However, in 
contrast to recently considered nonlinear PhC models, which assume transverse modulation 
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of the Kerr coefficient, we consider the case when only the refractive index is modulated, 
while the nonlinearity is uniform. Soliton dynamics in a ID version of this optical system, 
i.e., as a matter of fact, a multichannel nonlinear planar waveguide, was studied earlier in 
Ref. [17] (after the submission of the original version of this paper, a preprint has appeared 
which considers a very similar PhC model for the 2D case [18]). 

The present model is based on the rescaled Gross-Pitaevskii (GP) equation [19] for the 
wave function u(r,t), 



(in optical models, it is the NLS equation, t being the propagation distance). Here, V(r) is 
the external potential, which contains the parabolic-trap and OL terms, 



i.e., the OL period is normalized to be 7r (in the optical model, the parabolic trap is absent). 
The only dynamical invariant of Eq. (1) is the norm (number of atoms in BEC, or the 
beam's power in the optical model), N = J |w(r)| 2 dr. 

Analytical predictions for solitons can be obtained by means of the variational approxi- 
mation (VA), which was developed in nonlinear optics (see a review [20]), and successfully 
applied to BECs too [3,6,21], including BECs trapped in OL [22] (VA in the ID version 
of this model was earlier elaborated for the above-mentioned multichannel nonlinear opti- 
cal waveguide [17]). Following these works, we adopt the Gaussian ansatz for fundamental 
solitons, 



with an arbitrary negative constant /i and positive ones a and A (we set e > in Eq. (2), 
hence the central point of the ansatz must be chosen at the potential minimum, r = 0). In 
the case of BEC, /i is the chemical potential, and in optics —\i is the propagation constant. 
Substituting the ansatz into the Lagrangian L corresponding to Eq. (1), we perform spatial 



iu t + V 2 + V(r) + |u| 2 u = 



(1) 




(2) 




(3) 
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integration, and derive equations for a and A in the form dL/da = OL/dA — (// is not a 
variational degree of freedom, but rather an intrinsic parameter of the soliton family). 

Below, we display variational results for u = in the GP/NLS equation (1), which 
implies that the soliton has a size much smaller than that imposed by the parabolic trap. 
In the 2D case, the variational procedure amounts to an equation which determines a in 
terms of the norm, N = nA 2 /a [this expression corresponds to the ansatz (3)], and another 
equation that yields the chemical potential/propagation constant: 



Note that, if e — 0, Eqs. (4) yield iV = An. In the present notation, it is nothing else 
but a known variational prediction for the critical norm in the 2D NLS equation, which 
simultaneously is a universal value of the norm corresponding to the unstable 2D soliton. 
With e 7^ 0, Eq.(4) shows that the number of particles attains a minimum (threshold) value, 



at a = 1/2. This means that there is a threshold value of the norm which is necessary 
to create the 2D soliton; however, the threshold really exists only if it is positive. The 
expression (5) shows that the threshold exists for a relatively weak lattice, and it disappears 
if e exceeds the value e = e 2 /8 0. 92. 

It follows from Eq. (4) that the norm N cannot exceed the above-mentioned critical 
value, N cr = 4tt, which means that VA predicts a family of fundamental 2D solitons in the 
interval A^ thr < N < N CT (recall A^ hr = if e > e ). VA makes it also possible to predict 
stability of the solitons, on the basis of the Vakhitov-Kolokolov (VK) criterion [15], which 
applies if the solution to Eqs. (4) is obtained in a form /i = n(N): the necessary stability 
condition is dfi/dN < 0. 

A numerical solution of Eqs. (4), produces two different branches of the solution, see an 
example in Fig. la. An analytical consideration of the near-critical case, 1 — N/4ir — > +0, 
also yields two different solutions: 






(5) 
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fll (N)^-^ i _ N/ ^ +2ej, k(N) « {l-N/An)- In (2s)]- 1 . (6) 

Both the numerical and analytical solutions demonstrate that one branch \pi(N) in Eq. (6)] 
satisfies the VK criterion, and the other one [^(N) in Eq. (6)] does not. Direct simulations 
(see inset to Fig. la and more results below) confirm that the variational ansatz satisfying 
the VK criterion gives rise to a stable soliton, whose form is quite close to the predicted one. 

Qualitatively, the appearance of the family of stable solitons in the interval N thr < N < 
N cv in the presence of the lattice can be explained by the fact that it actually creates a 
nontrivial band (in terms of the corresponding values of /z), where the solitons may be stable 
(in the limit e — > the band degenerates into a single point). Note that appearance of a 
similar band explains the existence of stable 2D solitons in the OL-supported model with 
the positive scattering length [11]. 
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FIG. 1. (a) Numerical solution to Eqs. (4) for e = 0.4 and e = Eq = 0.92 (the latter value is 
the one at which AT thr vanishes). Inset displays evolution of the amplitude of a directly simulated 
solution initiated by the variational ansatz (3) with A and a taken as solutions of Eqs. (4), in the 
case of e = 0.92, N = 2ir, a = 1.3 (the VK-stable branch), (b) Dependence of the amplitude A 
of established 2D solitons vs. the initial norm N, as obtained from direct simulations of Eq. (1) 
with e = 0.4 and u = 0, starting with the configuration predicted by VA. The undulations in the 
dependence is a result of radiation loss in the course of the soliton formation. The dashed curve is 
the dependence A(N) as given by VA. 

The simple ansatz (3) does not take into regard distortion of the soliton's shape by the 
periodic potential; for this reason, VA is expected to be accurate if the soliton's size is not 
essentially larger than the lattice period, which is confirmed by comparison with numerical 
findings, see below. 

Direct simulations of Eq. (1) were performed by means of the multidimensional fast 
Fourier transform on grids of the size 256x256 or 128x128x128 in the 2D and 3D cases, 
respectively, with a time step At = 0.001. The domain size was \x,y\ < 4tt (in the figures 
we show only its central part). To emulate infinite boundary conditions (if the parabolic 
trap is switched off), absorbers were installed at borders of the integration domain, which 
eliminate linear waves emitted in the course of formation of solitons. 

A basic conclusion following from the simulations is that the initial waveforms, taken 
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as predicted by VA, always evolve into stable solitons (in both 2D and 3D cases). On the 
other hand, any input with the norm below the threshold predicted by VA (if e < e ), and 
some other input forms, taken very far from the variational prediction, decay into linear 
waves. This general result is illustrated by Fig. 1(b) from which we see that a delocalizing 
transition, manifested by the rapid drop of the amplitude, occurs at N cr /<Aii 0.62 for the 
case e = 0.4. By comparing this value with the minimum value attained by N/4tt along the 
VA existence curve of Fig. 1(a) at e — 0.4 (i.e. N min /4ir 0.58), we see that, although not 
perfect, the agreement is reasonable. Similar results are found also for other values of the 
amplitude of the optical lattice. We remark, however, that non-sharpness of the transition 
(see Fig. 1(b)) makes it difficult to determine the threshold with higher accuracy since, close 
to the threshold, long simulations are necessary to check whether a broad small-amplitude 
soliton is formed or not. Accurate study of the threshold norm versus e is presently under 
investigation and will be reported elsewhere. 

The presence of a lower threshold for the existence of the multi-dimensional solitons in 
the OL strongly resemble a similar phenomenon occurring in nonlinear lattices for intrinsic 
localized modes [23]. This analogy correlates with the fact that the above solitons can be 
strongly localized around one site of the lattice (see below). 

Actually, the lattice stabilizes the soliton not only against decay, but also partly against 
collapse. Indeed, in Fig. 1(b) the norm extends to slightly overcritical values (the exact 
value in 2D is N cx /An = 0.93 [15]) without provoking the collapse. An explanation for this 
comes from the fact that the chemical potential of the soliton lies in the gap of the excitation 
spectrum created by the lattice, thus enhancing their stability against the collapse or decay 
into linear waves. The mechanism for the formation of gap solitons in BECs with the OL 
and positive scattering length was identified as modulational instability of the Bloch states 
at the edges of the Brillouin zone [11,12]. 

Depending on the strength of the lattice potential and norm of the initial pulse, the 
emerging soliton occupies a single lattice cell (then we call it a "single-cell soliton"), or 
spreads itself over a few cells (a "multi-cell soliton"). Single-cell solitons always have a 
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larger norm than the multi-cell ones, and a border value of N, which separates the single- 
and multi-cell states, can be approximately identified. The relaxation of the 2D and 3D 
pulses into the stable shape proceeds with oscillations [see inset to Fig. 1(a)], with a higher 
frequency in stronger lattices. In fact, the input which is predicted by VA to be a stable 
soliton (according to the VK criterion) quickly evolves into a single-cell soliton whose shape 
is quite close to the predicted one, which is illustrated by inset to Fig. la. In contrast to 
this, the initial pulse taken as what is expected to be a VK-unstable soliton, according to 
VA, undergoes violent evolution, with large emission loss, which ends up by the formation 
of a multi-cell soliton. Full evolution pictures are not displayed here, as it is difficult to 
generate them in the 2D and 3D geometry without radial symmetry. 

Generic examples of the established single-cell and multi-cell 2D solitons, which were 
generated, respectively, by the VA-predicted VK-stable Gaussian ansatz in a moderately 
weak lattice potential, and the VK-unstable ansatz in a strong potential, are displayed in 
Fig. 2. A central peak of the soliton and satellites are clearly seen in the latter case. 




x x 
FIG. 2. (a) An established single-cell 2D soliton in a moderately strong lattice potential, found 

from direct simulations of Eq. (1) with e = 0.92. The initial configuration was taken as per the 

VA-predicted stable soliton, i.e., Eqs. (3) and (4) were used, with a = 1.3, the corresponding norm 

being N = 2tt. (b) An established multi-cell 2D soliton in a strong lattice potential, with e = 10. 

The initial configuration was taken as per the VA-predicted unstable soliton with a = 0.15 and 

N = 2tt [the same initial norm as in the case shown in (a)]. In the case (b), the norm drops to 

N = 3.18, i.e., w 50% of the initial norm is lost in the course of the formation of the stable soliton. 

The formation of solitons in the 3D lattice potential is, generally, similar to the 2D case, 
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although the relaxation time is shorter, due to the stronger interaction involving a larger 
number of adjacent cells. A typical example of a stable multi-cell 3D soliton, which forms 
itself in a strong lattice potential, is given in Fig. 3(a). 

Alongside the fundamental 2D soliton, in the absence of the lattice one can construct vor- 
tex solitons, of the form u = v(r) exp (— i(j,t + iS0), where r and are the polar coordinates, 
S is an integer vorticity ("spin"), and v(r) is a real function, which exponentially decays as 
r — > oo and vanishes ~ r s as r — > 0. As well as the fundamental (S = 0) soliton, the vortices 
are strongly unstable in the usual 2D NLS equation. Recently, stable vortices have been 
found in 2D (see reviews [24]) and 3D [25] models with non-Kerr optical nonlinearities. On 
the other hand, it was demonstrated [27] that stable vortices, solely with S = 1, exist in the 
discrete version of the usual cubic 2D NLS equation (which describes a bundle of nonlinear 
optical waveguides [26], or BEC trapped in a strong OL field [22]. Note that, unlike the 
isotropic NLS model, in ones with the axial symmetry broken by the lattice the vorticity is 
not a topological invariant, hence the very existence of such solutions is a nontrivial issue 
[27]. 

To generate vortex solitons, we simulated the 2D version of Eq. (1) with the initial 
condition u o (r,0) = Ar s exp(— or 2 ) exp(iS0), trying various values of A and a. Stable 
vortices with S > 2 could never be created (they quickly spread out); however, initial 
configurations with S = 1 readily self-trap into a stable vortex, which features a complex 
multi-cell organization, see Fig. 3(b). 
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FIG. 3. (a) An established 3D soliton formed in a strong lattice potential. The z = 
cross-section is shown for e = 10 and N = 10. (b) A typical example of a stable vortex soli- 
ton with S = 1 and N = 2tt, in the 2D lattice potential with e = 10. Inset shows the diagonal 
cross-section of the soliton; note that (as it should be) the field vanishes ~ r as r — > 0. 

All the above examples were obtained for uj = in Eq. (1), i.e., without the parabolic 
trap. Simulations with uj ^ have demonstrated that the external trap does not affect 
the 2D and 3D solitons in any conspicuous way, provided that the corresponding harmonic- 
oscillator length is essentially larger than the lattice period. 

The solitons found in this work are pinned by the lattice. In fact, a similar mechanism 
may stabilize solitons that partially keep the mobility of their free counterparts. To this 
end, one can use quasi- ID and quasi-2D lattices, in the 2D and 3D cases, respectively. This 
modification of the model will be reported elsewhere. 

Finally, we address experimental perspectives. First, creation of the initial waveform oc- 
cupying a single lattice cell can be done using the recently developed technique for patterned 
loading of BEC into optical lattices [28], which provides flexible control over placement of 
atoms in lattice sites. An initial waveform spread over multiple cells can be prepared by 
imposing an OL upon a condensate of a suitable size in the magnetic trap, with subsequent 
switching off the magnetic field. Changing the norm of the wave function, as supposed in the 
theory and numerical simulations, could be modeled by variation of the nonlinear coupling 
parameter (s-wave scattering length) via the Feshbach resonance. 

We note that the strengths of the OLs considered above are in the experimentally relevant 
range [10], e = 4- 20 in units of the recoil energy E rec = h 2 k 2 j /(2m), where k,L = 2n/X, 
and m and A denote the atomic mass and laser wavelength, respectively. We used the size 
of a unit cell d = A/2 = 0.425 //m, and the time 2m/(hk 2 L ) = 50 /is (for 85 Rb atoms and 
far detuned laser with A = 850 nm) as units of the length and time, respectively. Relevant 
experiments with can be performed with 7 Li or 85 Rb atoms featuring the negative scattering 
length in the ground state, which is amenable to large variations through the Feshbach 
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resonance. 

In conclusion, we have demonstrated that the GP/NLS equation with the attractive cubic 
nonlinearity and periodic potential, which describes BECs trapped in a 2D or 3D optical 
lattice, and (in the 2D case) an optical beam in the Kerr medium with a transverse periodic 
modulation of the refractive index, gives rise to stable solitons. In moderately weak and 
strong lattices, single-cell and multi-cell solitons were found, respectively, the former ones 
being accurately predicted by the variational approximation. A necessary condition for the 
formation of 2D solitons in a relatively weak lattice is that the initial norm of the field must 
exceed a threshold value. Stable 2D vortex solitons with 5 = 1 were found too. 
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